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I.  INTRODUCTION 


A  liquid  payload  in  a  spinning  projectile  can  cause  the  projectile  to  be 
unstable  for  certain  combinations  of  parameters  of  the  system,  even  though  the 
same  projectile  with  a  solid  payload  is  stable.  An  understanding  of  this 
phenomenon  requires  a  knowledge  of  the  wave  system  in  the  rotating  fluid.  The 
ability  to  support  waves  is  a  fundamental  property  of  a  rotating  fluid.  The 
angular  motion  of  the  projectile  excites  these  waves,  often  called  inertial 
waves.  In  this  study  a  numerical  experiment  generates  data  that  are  used  to 
validate  a  linear  theory  for  waves  in  a  rotating  fluid  contained  in  a 
cyl i nder. 

Only  linear  theories  are  available  to  analyze  the  effects  of  a  liquid 
payload  on  the  flight  of  a  spinning  projectile.  The  two  main  goals  of  such 
theories  are:  (1)  Calculate  the  moment  exerted  by  the  liquid  and  thereby 
determine  the  rate  of  growth  of  the  projectile  yaw-  (2)  Calculate  the 
frequencies  and  damping  of  the  predominant  modes  of  oscillation  of  the 
rotating  fluid.  For  a  cylinder  containing  a  liquid  in  solid-body  rotation,  an 
inviscid  theory  which  met  these  goals  was  presented  by  Stewartson,* 1  and  a 
significant  improvement  on  this  was  made  by  Wedemeyer,2 3  who  derived  viscous 
corrections.  Reference  3  presents  a  technique  for  achieving  (2)  for  a  liquid 
in  the  unsteady  spin-up  process  or  in  solid  body  rotation;  results  are  given 
there  only  for  the  latter  case. 

To  accomplish  (2),  the  first  step  is  to  linearize  the  Navier-Stokes 
equations.  A  small  perturbation  about  a  basic  state,  solid-body  rotation  or 
spin-up,  is  assumed,  and  these  perturbations  are  analyzed  into  modes.  The 
solution  to  the  eigenvalue  problem  thus  defined  yields  the  complex  eigenvalue* 
C  =  C^  +  iCj  where  C^  and  Cj  are  the  frequency  and  damping  of  the  wave, 

respectively.  With  a  modal  analysis,  the  no-slip  boundary  conditions  on  the 


* Definitions  of  terms  are  given  in  List  of  Symbols ,  p.  37. 

1.  K.  Stewartson ,  "On  the  Stability  of  a  Spinning  Top  Containing  Liquid ," 

J .  Fluid  Mech. ,  Vol.  5,  Part  4 ,  September  1959 ,  pp.  577-592. 

2.  E.  H.  Wedemeyer ,  "Viscous  Corrections  to  Stewartson' s  Stability  Criterion ," 
BEL  Report  No.  1325,  Aberdeen  Proving  Ground,  Maryland,  June  1966 

(AD  489687 ) .  Also  "Dynamics  of  Liquid-Filled  Shell:  Theory  of  Viscous 
Corrections  to  Stewartson' s  Stability  Problem ,"  BRL  Report  No.  1287, 

June  1965  (AD  472474). 

3.  C.  W.  Kitchens,  Jr.,  N.  Gerber,  and  R.  Sedney,  "Oscillations  of  a  Liquid 
in  a  Rotating  Cylinder :  Part  I.  Solid  Body  Rotation, "  Technical  Report 
ARBRL-TR- 02081 ,  June  1978  (AD  A057759). 
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end  walls  cannot  be  satisfied;  only  the  normal  component  of  velocity  vanishes 
there.  The  effect  of  using  the  "inviscid  type"  boundary  conditions  on  the  end 
wall  is  to  cause  an  error  in  CR  of,  at  most,  a  few  percent;  the  error  in  Cj 

can  be  more  than  a  factor  of  2.  Therefore,  a  correction  for  C  is  required. 

Such  a  correction  is  derived  in  Reference  3  for  perturbed  solid-body 
rotation.  Although  it  is  different  in  detail  from  that  derived  by  Uedemeyer,2 
the  necessary  analysis  proceeds  in  much  the  same  way.  The  attempts  in 
References  3  and  4  to  validate  the  linear  eigenvalue  theory  by  comparison  with 
experimental  data  were  not  definitive.  Validation  of  this  theory  is  the 
primary  objective  of  this  report. 

Recently,  a  finite  difference  technique  to  solve  the  Navier-Stokes  equa¬ 
tions  for  the  fluid  motion  in  a  rotating  cylinder  was  presented  by  Kitchens.5 
The  original  objective  of  the  work  was  to  study  the  axi symmetric  spin-up 
problem.  Since  this  technique  is  accurate  and  reliable,  it  provides  solutions 
where  the  accuracy  is  limited  only  by  the  grid  sizes.  Being  restricted  to 
axisymmetric  flow,  this  technique  can  only  yield  results  for  azimuthal  wave 
number  m  =  0.  For  the  projectile  case,  however,  m  =  1  is  the  wave  number  of 
interest.  We  believe  that  the  conclusions  we  reach  on  the  validity  of  the  linear 
theory  for  m  =  0  will  be  a  useful  guide  for  the  m  =  1  case.  Therefore,  we 
shall  use  the  technique  of  Reference  5  to  generate  data  with  which  to  validate 
the  linear  theory. 

Experimental  data  for  m  =  0  were  obtained  by  Aldridge;6  however,  a  sys¬ 
tematic  error  exists  in  the  reduction  of  the  data  to  frequency.  Recently, 
Aldridge  and  Stergiopoulos7  used  a  different  reduction  method  and  the  switch- 


4.  R.  D.  Whiting  and  N.  Gerber ,  "Dynamics  of  a  Liquid-Filled  Gyroscope : 

Update  of  Theory  and  Experiment , "  Technical  Report  ARBRL-TR-02221 , 

March  1980  (AD  A083886). 

5.  C.  W.  Kitchens ,  Jr.,  "Navier-Stokes  Solutions  for  Spin-Up  in  a  Fitted 
Cylinder, "  AIAA  Journal,  Vot.  18,  No.  8,  August  1980,  pp.  929-934.  Also 
Technical  RepormREWH-TR-02193,  September  1979  (AD  A077115). 

6.  X.  D.  Aldridge,  "Experimental  Verification  of  the  Inertial  Gscittations 
of  a  Fluid  in  a  Cylinder  During  Spin-Up , "  BRL  Contract  Report  No.  273, 
Aberdeen  Proving  Ground,  Maryland,  September  1975  (AD  A018797).  Also 
Geophys.  Astrophys.  Fluid  Dynamics,  Vot.  8,  1977,  pp.  279-301. 

7.  K.  D.  Aldridge  and  S.  Stergiopoulos ,  "Recovery  of  Complex  Eigen  frequencies 
of  Inertial  Waves  by  Linearized  Inversion  of  Viscous  Decay  Records , "  to 

be  submitted  for  publication.  See  also  S.  Stergiopoulos,  "An  Experi¬ 
mental  Study  of  Inertial  Waves  in  a  Fluid  Contained  in  a  Rotating 
Cylindrical  Cavity  During  Spin-Up  From  Rest,"  Ph.D  Thesis,  York  University, 
Toronto,  Ontario,  February  1982. 
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off  technique,  to  be  described  later,  and  obtained  frequencies  that  agree  with 
our  numerical  simulation  results.  In  these  experiments  and  our  finite 
difference  solutions  the  small  perturbation  is  a  sinusoidal  oscillation  of  the 
angular  velocity  of  the  cylinder.  Only  perturbed  solid-body  rotation  is 
reported  here;  spin-up  cases  have  been  done  and  will  be  reported  separately. 

Previous  experimental  determinations  of  the  frequencies  were  reported  by 
Fultz,8  using  an  oscillating  disk  on  the  axis  of  the  cylinder  to  excite  the 
waves  and  dye  to  observe  the  fluid  motion.  For  a  range  of  values  of  cylinder 
height  to  radius,  but  for  essentially  only  one  Reynolds  number,  he  found  that 
the  frequencies  agreed  closely  with  those  from  inviscid  theory.  The  theoreti¬ 
cal  study  by  Kudlick9  gave  viscous  corrections  to  the  inviscid  results  and 
expressions  for  the  wave  damping;  see  Section  V  for  further  discussion  of  this 
theory.  The  work  of  Kudlick9  and  Fultz8  is  discussed  by  Greenspan.10 

One  of  the  main  considerations  in  our  work  is  how  to  analyze  the  data  from 
the  finite  difference  solution.  Following  Aldridge,  we  use  pressure  differ¬ 
ences  on  the  axis  versus  time  as  our  diagnostic.  Of  the  doubly  infinite  set 
of  inertial  waves  generated  by  the  perturbation,  the  lower  modes  are  most 
pertinent.  If  the  damping  is  small  and  the  modes  are  well  separated,  we  can 
use  the  same  methods  for  analyzing  the  data  as  are  used  for  other  continuous 
systems,  which  are  based  on  simple  harmonic  oscillator  concepts.  For  several 
modes  these  conditions  were  not  satisfied,  and  for  these  cases  it  was 
necessary  to  use  Fourier  transform  and  digital  filter  methods.  For  Re  = 
0(102),  where  Re  =  Reynolds  number,  the  damping  in  the  system  was  so  large 
that  preprocessing  of  the  pressure  data  was  necessary. 

From  the  finite  difference  solutions  to  the  Navier-Stokes  equations,  the 
frequency  and  damping  were  obtained  for  several  modes  at  Re  =  102,  103, 
and  104,  mostly  for  aspect  ratio  unity.  For  the  two  higher  Re,  the  frequency 
agreed  with  to  within  +2%  for  the  k  =  2  mode  and  +3.5%  for  the  k  =  4  mode 
(mode  index  kKis  defined  in  (1)  and  (2)).  The  damping  agreed  with  Cj  to 

within  +5%  for  k  =  2  and  +6.5%  for  k  =  4.  It  appeared  that  closer  agreement 
with  Cj  could  have  been  obtained  if  the  accuracy  of  the  finite  difference 
solution  were  increased  by  decreasing  grid  sizes.  Since  the  correction  for  C 
is  of  a  boundary-layer  type,  its  validity  should  deteriorate  as  Re  decreases; 


8.  D.  Fultz ,  "A  Note  on  Overstability  and  the  Elastoid-Inertia  Oscillations 
of  Kelvin,  Solberg,  and  Bjerkness, "  J.  Meteorology,  Vol.  16,  April  1959, 
pp.  199-208. 

9.  M.  D.  Kudlick,  "On  the  Transient  Motion  of  a  Contained  Rotating  Fluid, " 
Ph.D.  Thesis,  Massachusetts  Institute  of  Technology,  1966. 

10.  H.  P.  Greenspan,  The  Theory  of  Rotating  Fluids,  Cambridge  University 
Press,  New  York,  NY,  1968,pp.  8lff. 
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at  the  lowest  Re,  larger  discrepancies  were  obtained.  The  conclusion  is  that 
the  linear  theory  with  correction  can  be  used  with  confidence  for  Re  as  low  as 
1,000. 


II.  THE  LINEAR  THEORY  WITH  CORRECTION 

A  brief  outline  of  the  theory  will  be  given,  mainly  to  introduce  notation 
and  concepts;  see  Reference  3  for  details.  The  unperturbed  container  spin 
rate  is  a  is  the  container  radius  and  a»c  is  its  half-height;  v  is  the 

liquid  kinematic  viscosity.  Length,  velocity,  pressure,  and  time  are  non- 

dimensional  ized  by  a,  aOQ,  pft2a2,  and  1/0Q,  respectively,  where  p  is  the 

liquid  density.  Dimensionless  nonrotating  cylindrical  coordinates  r,  e,  and  z 
are  used,  with  z  =  0  at  one  endwall;  dimensionless  time  is  denoted  by  t.  The 
coordinates  and  some  symbols  are  defined  in  Figure  1.  We  perturb  about  a 
basic  flow  of  solid-body  rotation  given  by 

U  =  0,  V  =  r,  W  =  0,  aP/3r  =  r, 

where  II,  V,  and  W  are  radial,  azimuthal,  and  axial  velocity  components, 
respectively,  and  P  is  pressure. 

★  ★  ★  ★ 

The  perturbation  pressure,  p,  and  velocities,  u,  v,  w  in  the  r,  e,  z 
directions,  satisfy  the  linearized  Navier-Stokes  equations.  At  this  stage 
there  is  no  restriction  on  Reynolds  number 

Re  =  a2/v 


unless  the  basic  flow  requires  one.  For  solid-body  rotation  there  is  no 
restriction.  We  assume  that  the  perturbation  can  be  represented  as  a 
superposition  of  modes: 


* 


u  =  u 

(r) 

COS 

Kz 

exp  [i  (Ct 

-  me)] 

(la) 

★  A 

V  =  V 

(r) 

cos 

Kz 

exp  [ 

i  (Ct 

-  me)] 

(lb) 

★  ~ 

w  =  w 

(r) 

sin 

Kz 

exp  [ 

i  (Ct 

-  me)] 

(lc) 

★  * 

p  =  p 

(r) 

cos 

Kz 

exp  [ 

i  (Ct 

-  me)] 

(Id) 

where 

K  =  k  tt/(  2c) .  (2) 


In  this  form  the  disturbances  are  complex,  the  real  parts  being  the  physical 

AAA  A 

quantities;  also  the  functions  u(r),  v(r),  w(r),  and  p(r)  are  complex.  The 
integers  m  =  0,  +1,  ...»  and  k  =  1,  2,  ...  are  azimuthal  and  axial  wave  num¬ 
bers,  respectively,  and  C  =  +  i  Cj  is  the  eigenvalue.  All  quantities  are 
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dimensionless.  The  natural  frequency  is  equal  to  CDft  and  the  decay  rate  of 

the  disturbance  is  equal  to  in  dimensional  form.  A  radial  wave  number, 

n,  is  determined  by  the  functions  of  r.  Since  the  numerical  simulation 
employed  here  requires  m  =  0  (axi symmetric  disturbance),  the  modes  can  be 
identified  by  the  pair  of  integers  (k,  n).  The  corresponding  values  of  Cp  and 
Cj  are  designated  with  the  subscripts  k,  n. 

The  functions  of  r  satisfy  a  6-th  order  system  (complex)  of  ordinary 
differential  equations  which,  together  with  homogeneous  boundary  conditions  at 
r  =  0  and  r  =  1,  define  the  eigenvalue  problem.  The  no-slip  boundary  condi¬ 
tions  are  satisfied  at  the  side  wall,  r  =  1.  Since  the  equations  have  a 
singularity  at  r  =  0,  the  form  of  the  solution  is  determined  analytically  near 
r  =  0;  the  boundary  conditions  there  must  be  derived  on  the  basis  of  conti¬ 
nuity  and  single-valuedness.  For  Re  >  0(103)  special  care  is  required  in  the 
numerical  integration.  The  process  of  orthonormalization  is  used  to  cope  with 
the  phenomena  of  growing  solutions  and  loss  of  linear  independence. 

<n> 

The  boundary  conditions  at  the  end  walls  (z  =  0,  2c)  are  u  =  v  =  w  =  0. 

* 

The  modal  decomposition  form  assumed  in  (1)  allows  only  the  condition  w  =  0  to 

be  satisfied;  the  no-slip  condition  on  the  end  walls,  u  =  v  =  0,  is  not 
satisfied.  To  account  for  this  a  boundary- layer  type  of  correction  is  made  to 
the  solution;  it  is  described  in  Reference  3. 

The  correction  is  essentially  a  solution  to  the  perturbation  boundary 
layer  equations.  The  boundary  layer  is  introduced  to  correct  for  the  neglect 
of  the  no-slip  boundary  conditions  at  the  end  walls.  From  the  solution  of 
these  boundary  layer  equations  a  complex  displacement  thickness,  6c,  is 
determined.  The  K  of  (2)  is  replaced  by  K£  =  (ku/2)/(c  -  6c).  The  eigenvalue 

problem  is  solved  again  to  yield  the  corrected  C.  The  solution  to  the 

perturbation  boundary- layer  equations  is  correct  to  0(Re”^/^);  therefore,  a 
restriction  to  "large"  Re  now  exists. 

Two  questions  can  be  raised:  (1)  Is  the  corrected  C  adequate  for  present 
application.  (2)  What  is  the  lower  limit  on  Re  for  a  given  accuracy  of  C? 

Ordi nari ly,  one  would  seek  experimental  results  to  answer  these  questions. 
However,  there  was  a  scarcity  of  experimental  data  at  the  time  of  the  early 
computations  of  C.  The  available  measurements  of  Aldridge6  consisted  only  of 
frequency  mainly  for  spin-up,  and  the  range  of  parameters  in  his  experi¬ 
ments  was  not  sufficient  for  our  purposes.  Therefore,  we  resorted  to  a  numeri¬ 
cal  simulation  to  obtain  frequency  and  damping  with  which  C  could  be  compared. 


III.  FINITE  DIFFERENCE  SOLUTION  TO  NAVIER-ST0KES  EQUATIONS 
FOR  A  ROTATING  CYLINDER 

Kitchens5  developed  a  finite  difference  technique  for  solving  the  Navier- 
Stokes  equations  for  axisymmetric  flow  in  a  rotating  cylinder.  This  solution 
was  set  up  especially  for  the  spin-up  flow  but  can  be  adapted  to  other  flows 
by  changing  the  boundary  and/or  initial  conditions.  Because  of  the  special 
problem  considered  in  Reference  5,  symmetry  about  the  midplane,  z  =  c,  was 
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assumed;  this  restriction  can  easily  be  removed  but  will  be  retained  here; 
thus,  only  even  k  modes  are  possible. 

Coordinate  transformations  allow  concentration  of  grid  points  near  the 
boundaries.  A  stability  criterion  giving  the  allowable  step  size.  At,  is 
tested  at  each  step  of  the  calculation.  Seven  other  parameters  determine  the 
convergence  and  accuracy  of  the  solution:  the  numbers  of  points  in  the  r  and 
z  directions,  relaxation  constants  for  each  of  the  three  dependent  variables, 
and  two  coordinate  stretching  constants.  This  method  has  proved  to  be 
accurate  and  reliable.  For  spin-up  from  rest,  convergence  was  obtained  for  Re 
up  to  5  x  104  with  moderate  computation  time;  convergence  was  not  achieved  for 

Re  =  in5.  For  the  perturbation  problem  calculations  for  larger  Re  have  been 
made. 

Since  we  wished  to  validate  linear  theory,  calculations  were  made  with  a 
small  perturbation  on  solid  body  rotation.  The  time  history  of  angular 
velocity  prescribed  for  the  cylinder  is 

n(t)/nQ  =  1 

=  1  +  esin  cot 
=  1 

For  a  range  of  values  of  e  near  e  =  0.05  it  was  verified  that  the  solution 
depended  linearly  on  e.  The  forcing  function  ft(t)  induces  a  wave-type 
response  in  the  fluid.  The  velocity  components,  pressure,  and  stream 
function,  which  are  computed  at  all  grid  points,  are  available  to  analyze  the 
response.  A  choice  must  be  made  as  to  the  flow  variable  and  the  particular 
aspect  of  its  variation  to  be  used  in  extracting  the  wave  parameters. 

A  physically  significant  diagnostic  is  the  pressure  difference  on  the  axis 
of  the  cylinder.  A  derivation  of  the  formula  for  pressure,  p(r,z,t)  is  given 
in  the  Appendix.  In  the  experiments  of  Aldridge  two  pressure  differences 

APj  =  p  (n,  c,  t)  -  p  (o,  n,  t) 


t  <  o 
0  <  t  <  t 

( 

t  >  t 


and 


UP\/2  =  P(°»  c/2 ,  t)  -  p  (0,  0,  t) 


were  used  to  extract  data  for  k  =  2  and  k  =  4,  respectively.  In  the  reduction 
of  our  data  it  was  found  that  Ap^  did  not  clearly  separate  k  =  2  and  k  =  4  . 

However, 

Apa  =  0.5  APi  -  Ap1/? 

effectively  separates  k  =  4.  The  reason  for  this  can  be  seen  if  it  is  assumed 
that  the  pressure  at  r  =  0  can  be  expanded  in  a  series  of  eigenfunctions 
obtained  from  (1): 
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(3) 


p(0,  z,  t)  =  E  A.  cos  Kz  exp  (iC.t) 

L,  _  O  A  ^  fl 

^  H-j  ••• 

n  =  1,  2,  ... 

where  A^n  are  constant  coefficients.  Front  this 

APi  =  -2  E  A2n  exp  (iC2nt)  +  ...  (4) 


Ap 


1/2 


=  -  E 

n 


AZn  ex'’ 


<1C2nt) 


E 

n 


4n 


exp  (iC4nt)  + 


(5) 


Apg  -  3  A^n  exp  ( iC4r|t)  +  ...  (6) 

are  obtained,  where  the  omitted  terms  start  with  k  =  6.  It  is  clear  from 
these  expressions  why  Ap^2  does  not  separate  out  the  k  =  4  mode  whereas 

Apa  does;  this  is  most  evident  if  C2i  =  C41. 

An  example  of  the  time  history  of  Api  is  shown  in  Figure  2.  The  analysis 
of  such  records,  which  yields  the  frequency  and  damping  of  the  inertial  waves, 
is  discussed  next. 


IV.  ANALYSIS  OF  THE  DATA 


The  response,  Ap^  vs  t,  illustrated  in  Figure  2  is  simple  enough  so  that 

ideas  from  simple  harmonic  oscillator  analysis  can  be  used;  for  this  case  e  = 
0.05,  a)  =  1.25,  t  =  70.4.  The  transient  response,  caused  by  initiating  the 

perturbation,  lasts  until  about  t  *  50  and  the  steady-state  response  exists 
for  50  <  t  <  t  =  70.4.  For  t  >  t  the  free  response  decays  approximately 

exponentially.  The  step  size  in  this  example  is  At  =  0.125  so  that  there  are 
forty  or  more  points  per  cycle  of  response;  in  all  cases  this  number  was  never 
less  than  twenty. 


Amplitude  response  curves,  obtained  by  plotting  max  |api|  and  max  | Ap  |  of 

the  steady  state  vs  w,  are  shown  in  Figure  3  for  Re  =  1.0  x  103  and 
1.0  x  10\  and  c/a  =  1.0.  The  nondimensional  flow  solution  is  completely 
characterized  by  Re  and  c/a.  The  values  of  u>  at  the  peaks  on  the  curves, 


a)^,  approximate  the  frequencies  of  inertial  oscillations,  f jcn;  the  values  of  n 
increase  from  n  =  1  as  to  decreases.  From  simple  harmonic  oscillator  analysis 
the  error  in  determining  fkn  from  the  peak  location  is 
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(f.  -  0)  )/f.  «•  (a  /f.  )2 

kn  m  kn  kit  kn 


where  is  the  damping  of  the  (k,  n)  mode.  The  sharp  peaks  have  small 

values  of  a^/f^.  F°r  =  1»0  x  103  the  second  peak  for  k  =  2  is  not  well 

defined,  and  the  k  =  4  curve  has  no  peaks;  f  cannot  be  determined  from  the 
amplitude  response  curves  in  these  cases. 


For  the  type  of  time  response  illustrated  in  Figure  2, the  damping, 
a2i,  can  be  determined  from  the  logarithmic  decrement  of  the  peaks  of  the  free 
response  for  t  >  tc»  If,  however,  two  or  more  modes  are  mixed  in  the  free 

response,  the  above  method  cannot  by  itself  be  used  to  determine  a.  The 
presence  of  two  mixed  modes  can  be  surmised  from  the  pressure  histories;  a 
definite  indication  can  be  obtained  by  plotting  the  instantaneous  stream¬ 
lines.  Since  this  is  a  more  time-consuming  process,  it  is  not  done  routinely. 

4 

These  are  shown  in  Figure  4  fo£  Re  =  10  ,  c/a  =  1,  e  =  0.05  at  four  values 
of  a).  In  Figures  4a  and  4c  co  =  u>  in  the  amplitude  response  curve  for  k  =  2 

in  Figure  3,  and  well-defined  cells  are  plotted  for  n  =  3  and  n  =  2,  respec¬ 
tively,  showing  the  presence  of  single  modes.  In  Figure  4b  oj  =  0.7,  a  value 
between  two  peaks  of  the  amplitude  response  curve.  As  expected,  well-defined 
cells  are  not  obtained,  showing  the  presence  of  mixed  modes.  In  Figure  4d  co  = 
1.05,  which  is  near  a  minimum  in  the  k  =  2  and  a  maximum  in  the  k  =  4 
amplitude  response  curves,  and  well-defined  cells  are  obtained  for  the 
k  =  4,  n  =  3  mode. 


For  Re  =  10 3  a  Apj  vs  t  response  with  two  modes  present  in  the  free 

response  is  given  in  Figure  5.  Only  the  response  after  cutoff,  t  >  t  = 

67.875,  is  used  in  the  analysis.  The  parameters  other  than  t  for  Figure  5 

are  the  same  as  those  for  Figure  2  except  that  u>  =  0.88,  which  is  the  value 
of  at  the  broad  maximum  in  the  Apx  curve  in  Figure  3  for  Re  =  103.  The 

lack  of  a  sharp  peak  indicates  that  the  elementary  method  may  fail. 

Figure  5  suggests  that  there  are  two  or  more  modes  mixed  in  the  transient 
response  data,  and  this  is  verified  by  the  Fourier  transform  of  this  data 
shown  in  Figure  6.  The  magnitude  of  the  transform,  |F|,  is  plotted  vs  the 

frequency,  v^,  in  Hz.  The  dimension  of  assumes  that  J2  =  1  rad/s.  A  fast 

Fourier  transform  computer  program  is  used  to  calculate  |F|.  In  obtaining  |F 
for  this  case,  data  from  67.875  <  t  <  125.0  are  used,  giving  458  points.  The 
proqram  requires  a  number  of  points  equal  to  a  power  of  2;  to  obtain  a  con¬ 
venient  interval  Avf  for  calculating  |F|,  zeros  are  added  to  the  sequence  of 

data  points  to  give  4,096  points.  The  transform  is  shown  in  Figure  6  for 
vf  <  0.4  Hz,  although  it  was  calculated  for  0  <  Vf  <  4.0  Hz;  the  part  beyond 

0.4  Hz  is  of  no  interest. 


Aside  from  the  solution  to  the  Navier-Stokes  equations,  the  Fourier 
transform  depends  on  the  number  of  data  points  and  their  location  in  the 


14 


complete  data  set.  Deleting  points  from  the  end  of  the  interval  has  little 
effect  on  ]F|,  but  deleting  points  at  the  beginning  may  have  a  noticeable 
effect.  If  the  data  starts  at  t  =  tc,  it  may  contain  enough  of  the  steady- 

state  response  so  that  a)  appears  in  the  spectrum  or  influences  it.  If  the 
initial  data  point  is  several  time  units  greater  that  tc  and  there  is  a  mode 

with  relatively  large  damping,  that  mode  may  not  appear  in  the  spectrum  or  it 
may  be  distorted.  In  our  treatment  the  optimum  choice  of  initial  point  was 
determined  empirically.  If  the  initial  point  for  the  example  in  Figure  6  is 
taken  to  be  t  =  75.5,  the  general  shape  of  the  curve  is  the  sane  and  the  two 
frequencies  in  the  spectrum  are  unchanged,  but  |F|  is  reduced  by  a  factor  of 
3,  and  some  additional  small  peaks  appear. 

To  obtain  an  accurate  value  of  the  higher  frequency  in  the  spectrum,  an 
a)  closer  to  f21  than  to  f 22  should  be  used  to  generate  the  data;  an  initial 

guess  can  be  obtained  from  inviscid  theory.  Then  f 22  is  obtained  using  an 

a)  closer  to  f  22.  For  the  particular  case  shown  in  Figure  6,  f 2\  and  f 22  are 

both  determined  to  within  ±\%  by  the  spectrum  for  m  =  0.88;  but  this  is  not 
always  the  case.  In  general,  the  accurate  evaluation  of  the  eigenf requencies 
requires  that  they  be  determined  in  decreasing  order. 

The  damping  is  obtained  after  the  data  is  processed  by  a  digital  filter. 
Usually  a  band  pass  filter  is  used  with  the  high  and  low  band  pass  limits 
determined  by  a  frequency  band  around  a  peak  in  the  |F|  curve.  The  result  of 
applying  a  band  pass  filter  with  limits  (0.11,  0.15)  Hz  to  the  data  in  Figure 
5  is  shown  in  Figure  7.  The  damping  a22  is  obtained  from  the  logarithmic 

decay  of  the  maxima  and  minima  in  this  filtered  data. 

The  filter,  however,  introduces  distortion  at  the  beginning  and  end  of  the 
data  range.  In  some  cases  only  a  few  maxima  and  minima  are  free  from  distor¬ 
tion  so  that  spurious  values  of  are  obtained  if  all  are  used.  Discarding 

all  but  the  few  not  distorted  increases  the  error  in  o^.  This  error  is  a 

function  of  Re,  k,  and  n.  The  number  of  undistorted  maxima  and  minima  could 
be  increased  by  judiciously  adjusting  the  parameters  of  the  system. 

Another  difficulty  occurred  when  the  above  process  was  applied  to  Re  =  100 
and  200  cases.  The  damping  was  so  large  that  accurate  values  for  f21  could 

not  be  obtained  and  a  second  peak  in  |F|,  to  determine  f 22,  could  not  be 
found.  A  method  of  preprocessing  the  data  was  developed  to  overcome  this. 

It  will  be  illustrated  for  k  =  2,  n  =  1,  2. 

It  is  assumed  that  the  data  can  be  approximated  by  the  form  of  expansion 
given  in  (4)  plus  a  contribution  from  the  noise  in  the  data,  n(t),  arising 
from  the  finite  difference  approximation.  Therefore,  with  obvious  changes  in 
notation  in  (4),  the  Ap  ^  data  can  be  expressed  as 


APi  =  a1e”0t2ltsin  f 


21t  +  a2e~a22tsin  f 


22 


t  + 


+  n(t) 
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where  t  =  t  -  t£  and  ai  are  related  to  Ai  and  the  other  constants.  Let  a21  be 
a  first  approximation  to  the  unknown  a2i  such  that  a21  <  a21;  if  the 
inequality  is  not  satisfied,  this  will  be  apparent  in  the  preprocessed  data. 

The  APi  data  is  modified  by  multiplying  it  by  ea2lt.  Thus, a  new  time 
series  is  generated,  giving  the  preprocessed  data: 


>a2 


APi  =  aie'(a21  -  «2l)tsin  f2]i  + 


a  2e 


" (  a22  ”  a21 


'sin  f22t  +  ...  +  ea2lt  n(t) 


(7) 


The  noise  is  negligible  except  near  the  end  of  the  record  where  Ap:  is  small. 

Since  n(t)  is  not  an  exponentially  damped  function,  the  last  term  in  (7)  is 

amplified  and  diverges  for  increasing  t.  However,  this  does  not  prevent  the 

determination  of  f 21  and  f 22  from  the  Fourier  transform  of  the  preprocessed 

data  because  the  eigenf requencies  are  less  than  1  Hz,  whereas  the  noise 
contributes  to  higher  frequencies,  typically  2  Hz  and  greater.  From  (7),  the 
effective  damping  is  reduced,  enabl i ng  f21  and  f22  to  be  determined  if 

“21  <  a2i*  The  data  will  diverge  for  t  >  0  if  a21  >  a21.  In  that 

case  a21  must  be  decreased  and  the  process  repeated.  It  is  feasible  to  try 

several  values  of  a2i  since  obtaining  the  preprocessed  data  is  simple  and 

inexpensive.  The  preprocessed  data  are  treated  as  before  to  yield  a21  - 

a2i  and  a22  -  a2i« 

An  example  of  the  result  of  preprocessing  is  shown  in  Figure  8.  The 
Apx  data,  for  Re  =  100,  are  in  the  form  of  a  highly  damped  sine  wave.  The 

Fourier  transform  of  this  data,  shown  in  Figure  8a,  has  only  one  peak  at  = 

0.177  Hz,  corresponding  to  f21.  The  transform  of  the  preprocessed  data 

with  ct21  =  0.30,  shown  in  Figure  8b,  has  two  peaks.  The  peak  at  the  higher 

frequency  has  shifted  to  v^-  =  0.189  Hz,  which  is  a  more  accurate  value 

for  1 21  than  0.177  Hz.  The  peak  at  the  lower  frequency  occurs  at  =  0.089 

Hz,  giving  f 22.  Values  for  f 2^  and  f 22  to  within  3%  of  the  above  values  are 

obtained  with  a21  =  0.36.  These  two  values  of  a21  can  be  conpared  with  a2i  = 

0.40,  the  value  finally  determined.  Preprocessing  thus  enables  the  f. 

k  n 

and  a)<n  to  be  determined  for  Re  =  100  and  200;  however,  their  accuracy  is  less 
than  that  for  higher  Re. 
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V.  RESULTS 


The  fkn  and  were  determined  by  the  procedures  described  in  Section  IV 

for  102  <  Re  <  5  x  104,  c/a  =  1,  n  =  1,  2,  3,  and  k  =  2,  4.  Some  points  in 
this  parameter  space  are  not  included  here.  These  results  are  compared  with 
the  and  Cj  obtained  from  the  linear  theory. 


In  Figure  9  the  f0  are  compared  with  the  ei genf requenci es  CD  .  For  Re 
1  R2n 

=  ID3  and  n  =  1,  2  the  f„  and  Cn  are  both  essentially  the  same  as  the 

R2n 

values  calculated  for  inviscid  flow;  however,  significant  difference  is  found 
for  f  23  at  Re  =  10 3. 


For  Re  <  10 3  the  differences  between  f^n  and  increase.  The  results 

for  Re  =  100,  200  are  shown  by  bars  indicating  the  variation  obtained  as  the 
various  parameters  of  the  analysis  process  are  changed;  for  the  other  symbols 
the  variation  is  +1%.  The  f 22  appear  to  be  nonmonotonic  at  the  lower  Re. 

The  Cj?  are  nonmonotonic  to  a  slight  degree.  In  the  application  of  the 

calculated  frequencies  to  the  liquid-filled  projectile,  an  accuracy  of  a  few 
per  cent  is  required. 


In  Figure  10  the  01^  are  compared  with  .  They  agree  to  within  ±5%  for 
Re  >  10 3,  but  the  deviations  are  considerably  more  for  lower  Re,  the  largest 


being  for  a  at  Re  =  100;  a  is  definitely  nonmonotonic.  For  n  =  1  the 
22  22 

result  from  the  theory  of  Kudlick9  is  also  shown.  This  theory  uses  a  matched- 

-1/2 

asymptotic  expansion  in  powers  of  Re  ',  but  only  the  lowest  term  in 
-1/2 

Re  is  calculated.  The  straight  line  from  Kudlick 's  theory  fairs  into  the 

curve  for  CT  at  about  Re  =  105,  but  it  deviates  increasingly  for  Re  <  105. 

1  21 

If  the  curve  of  C.  vs  Re  for  n  =  1,  from  the  linear  theory,  is  fitted  to 

1  -1/2  -1  04 

powers  of  Re,  starting  from  Re  ,  the  next  term  is  found  to  be  Re  *  , 


indicating  that  the  next  term  in  Kudlick's  expansion  is  necessary  to  obtain  a 
result  which  agrees  with  those  from  the  numerical  simulation.  The  presence  of 
log  Re  in  the  expansion  might  not  be  detected  from  this  curve  fit. 


The  same  type  of  results,  for  frequency  and  damping  vs  Re,  are  shown  in 

Figures  11  and  12,  respectively,  for  k  =  4.  For  Re  =  10  3’  f4n  -  are 

n4n 

slightly  larger  than  for  the  k  =  2  case  whereas  cu  -  C.  are  negligible,  as 

n  *4n 

for  k  =  2.  The  differences  between  the  results  of  the  Kudlick  theory  and 
those  of  either  the  linear  eigenvalue  theory  or  the  numerical  simulation  are 
greater  for  the  k  =  4  case  than  those  shown  in  Figure  10  for  k  =  2. 
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VI.  CONCLUSIONS 


By  means  of  a  numerical  simul ation, the  results  of  a  linear  eigenvalue 
theory  with  correction  have  been  validated.  Although  no  results  without  the 
correction  were  discussed,  the  calculations  showed  that  errors  of  several 
percent  in  CR  and  up  to  a  factor  of  2  in  Cj  are  obtained  if  the  correction  is 

not  applied.  With  correction,  the  predicted  C^  agree  with  the  f  to  within  +2% 

for  the  k  =  2  mode  and  3.5%  for  the  k  =  4  node;  the  CT  agree  with  the  a  to 

1  3 

within  +2%  for  the  k  =  2  mode  and  6.5%  for  the  k  =  4  node,  for  Re  >  10  .  For 

Re  <  103  extraction  of  f  and  a  from  the  pressure  data,  using  Fourier  transform 
and  digital  filter  techniques,  becomes  more  difficult,  and  the  uncertainty  in 
the  results  increases. 
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Figure  1.  Cylinder,  Coordinates,  and  Notation 


Figure  2.  Typical  Pressure  History  -  Ap^  vs  t  for 

Perturbed  Solid  Body  Rotation  with  Single 
Mode  Excitation  (Re  =  1.0  x  103,  c/a  =  1.0, 
s  =  0.05,  a)  =  1.15,  t  =  70.375). 
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Figure  3.  Amplitude  Response  Curve  -  max  j Ap |  vs 
Forcing  Frequency  (c/a  =  1.0,  e  =  0.05). 
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a)  =  0,6 

CR  (k=2,  n=3)  =  0.585 


Z 


to  =  0.7 


Figure  4.  Instantaneous  Meridional  Streamlines  for  Four  Perturbation 
Frequencies  (Re  =  lxlO4,  c/a  =  1.0,  e  =  0.05). 
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u  =  0.8 


CR  (k=2,  n  =2)  =  0.812 


Z 


(d) 


r 


a)  =  1.05 

CR  (k=4,  n=3)  =  1.05 


Figure  4.  Instantaneous  Meridional  Streamlines  for  Four  Perturbation 
Frequencies  (Re  =  lxlO4,  c/a  =  1.0,  e  =  0.05)  (continued). 
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Figure  5.  Pressure 

With  Two 
1.0,  e  = 


History,  Ap^  After 

Dominant  Modes  (Re 
0.05,  co  =  0.88,  tc 


Cut-Off  for  a  Case 

=  1.0  x  10^,  c/a  - 
=  67.875). 


vf  (Hz) 


Figure  6.  Fourier  Transform  of  the  Pressure  History 
Shown  in  Figure  5  (Re  =  1.0  x  10^,  c/a  = 
1.0,  e  =  0.05). 
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Figure  7.  Band  Pass  Filter  of  the  Pressure  History 

O 

Shown  in  Figure  5  (Re  =  1.0  x  10  ,  c/a  = 
1.0,  e  =  0.05,  t  =  67.875,  0.11  <  v.  S 
0.15).  c  r 


24 


o.oto 


|F| 


0  009 


0.009 

0.004 

0.002 

0.000 


■j 

- 

\ 

\ 

V 

- 

\ 

- r  ii  i 

[ti-vt 

pm  i  1 

i  t  i  r 

l'1'  I1 1  TT  ‘1  ""j 

1  1  1  1  1 

0.00  0.03  0.10  0.13  0.20  0.23 


(a) 


vf  (Hz) 


Figure  8.  Comparison  of  Fourier  Transforms  (a)  Before 
and  (b)  After  Preprocessing  the  Data  (Re  = 
100.0,  c/a  =  1.0,  e  =  0.05,  oj  =  0.83,  t  = 

40.0,  40.0  <  t  <  53.0,  a  =  0.3). 
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Figure  9.  Comparison  of  Eigenfrequencies  From  Linear  Theory  With  Correction,  CR,  and 
Those  From  Numerical  Simulation,  f  (c/a  =  1.0,  k  =  2). 
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Figure  10.  Comparison  of  Wave  Dampings  From  Linear  Theory,  Cj,  With  Those  From 
Simulation,  a,  and  From  Kualick  Theory  (c/a  =  1.0,  k  =  2). 


Figure  11. 


Comparison  of  Eigenfrequencies  From  Linear  Theory  With 
Correction,  Cp,  and  Those  From  Numerical  Simulation, 

f  (c/a  =  1.0,  k  =  4). 
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1.00 


F-qure  12.  Comparison  of  Wave  Dampings  From  Linear  Theory,  Cj,  With 

Those  From  Numerical  Simulation,  a,  and  From  Kudlick 
Theory  (c/a  =  1.0,  k  =  4). 
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APPENDIX:  DETERMINATION  OF  PRESSURE  FROM  STREAM  FUNCTION 


The  nondimensional  Navier-Stokes  axial  momentum  equation  (e.  g., 
(3.33c)  of  Reference  A-l)  for  incompressible  axi symmetric  flow  is  integrated 
with  respect  to  z  at  r  =  0: 


p(0,z,t)  -  p(n,nft)  =  C  i( 1/Re)  [wpr  +  (1/r)  wp  +  - 

Ju 


( A.  1 ) 


Cwt  +  vvwz]|dz, 

where  subscripts  denote  partial  derivatives.  Kitchens5  did  not  calculate  w, 
but  rather  the  stream  function  ij;.  Ry  definition 


w  =  -  i|>p/r.  (A.?) 

This  shows  how  p  is  obtained,  in  principle.  Its  evaluation  from  the 
finite  difference  solution  is  outlined  next.  It  is  inferred  from  the  equa¬ 
tions  and  boundary  conditions  in  Reference  S  that  t| >  is  an  even  function  of  r 
and  \p  (r  =  n)  =  0.  Then  ip  is  expanded  in  the  following  form  near  the  axis: 


♦  (r,z,t)  =  a,,  (z,t)r  +  a.  (z,t)r  + 


(A. 3) 


It  follows  from  (A.?)  that 

w  (r,z,t)  =  -2  a2  (z,t)  -  4  a4  (z,t)r?  +  ...  (A. 4) 

Evaluating  the  terms  on  the  right-hand  side  of  (A.l)  at  r  =  0  by  means  of 
(A. 4) ,  we  obtain 

9  rz 

p(n,z,t)  =  p(n,n,t)  =  -2.(a?(z,t)V  +  ?.  “  l  a£  U,t)d£  - 

•Ml 

(16/Re)  I  a4  U,t)d£  -  (A. 5) 

Jo 

(2/Re)  [(9a?/9z)z  _  z  -  (aa2/az)z  =  ()'J. 


A-l.  H.  Schlichting ,  Boundary  Layer  Theory ,  McGraw-Hill  Book  Co. ,  New  York , 
NY,  1960.  ~~  ^ 


On 

Ov> 


The  functions  a^z, t)  and  a^(z,t)  are  obtained  from  the  finite-difference 

solution5  at  those  values  of  z  at  which  vp(r,z,t)  is  computed.  A  grid  point  in 

the  finite-difference  solution  is  denoted  by  r^,  z^,  and  the  corresponding 

stream  function  by  j  =  1  for  r  =  n,  and  I  =  1  for  z  =  0.  The 

coefficients  an<^  a4»  evaluated  at  z  =  z^,  are  denoted  by  anc1  a4^> 

respectively.  If  we  retain  only  two  terms  in  the  expansion  (A. 3),  then  in 
terms  of  p?  ^  av<j> 


a?Z~^r3/r2^  ^2Z  "  ^ r?/r l)  ^3£~^r3  r2 


(A. 6) 


2  2- 


4£ 


—  ( 1  /  r ^ ^  ( 1/  r3 )  ^3 ^ T  r 3  r^  1 


Although  the  a^  and  a ,,  are  evaluated  only  at  z^,  they  are  continuous 
functions  of  z,  and  the  integrals  in  (A.b)  are  approximated  by  numerical  inte¬ 
gration  employing  the  values  of  a^  and  a^  .  Trapezoidal  integration  is  used 
for  z  =  Z£  and  Simpson  integratior  is  employed  for  z  >  z?.  For  a?,  e.g.,  the 
integral  is 


L>  ■ ! 


a?  (dz/dn)  dn. 


0 


where  n(z)  and  dz/dn  are  obtained  from  the  stretching  transformation  of  (17) 
in  Reference  5.  Then 


!t  1  >*-2  +  (4n/6)  r(a2>t-2  (dz/d">*-2  *  4  (a2>t-l  <dz/d,'>«-l  * 


a2*  (dz/dr'>P> 


where  An  is  the  constant  interval. 

The  3a,,/ 3z  terms  are  obtained  from  the  expression 


(A. 7) 


a2  (z)  =  a2 Jt  +  C1  <Z'ZP  *  C2  (z-zt)  • 


Then 

3a 


2)  _  r.<a?)3  -  (yzp  -  rta2)b  -  (a2)tl  <za-zt)? 

82  4- z,'  1  "  (za  -  V  (zb  -  zt>  (zb  -  za> 


(A.R) 


(A. 9) 
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where  z,  and 

d 

Finally,  the 
approximated 


where  At  is 


are  grid  points  either  adjacent  to  or  one  point  removed 
time  derivative  of  a  function,  g,  beginning  at  t  =  ?At  is 

by 

(dy/dt)t=t_  =  [3  gi  -  4  gi_1  +  gi_23/(2At) ,  (A 

the  constant  time  interval  of  the  calculation. 


10) 
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LIST  OF  SYMBOLS 


cross-sectional  radius  of  cylinder 

coefficient  of  kn  tern  in  z,t  expansion  of  p(0,z,t)  in  (3) 
half-height  of  cylinder/a 

=  CR  +  i  Cj,  complex  ei genf requency  of  rotating  fluid 

decay  rate/QQ  of  inertial  oscillation  computed  by  linear  theory 
frequency/ of  inertial  oscillation  computed  by  linear  theory 

frequency/ of  inertial  oscillation  determined  by  numerical 
simulation 

magnitude  of  Fourier  transform  of  pressure  history  data 

index  in  r-di recti  on  of  grid  point  in  finite-difference  solution 

axial  wave  number  ((la)  -  (Id)) 

=  kir/(2c) 

=  ( k-n/2 )  /  (c- 6) 

index  in  z-direction  of  grid  point  in  finite-difference  solution 

azimuthal  wave  number  (see  (la)  -  (Id)) 

radial  wave  number 

=  P  +  p,  total  pressure/ ( pa 2a02) 

perturbation  pressure/( pa2fy,2) 

radial  variation  of  p  (see  (Id)) 

unperturbed  pressure/( pa2aQ2) 

radial  coordinate 

=  ^  a2/v,  Reynolds  number 

time  x 

cutoff  time  x  of  perturbation  on  angular  velocity  of  cylinder 


U,  V,  w 

*  ★  ★ 
U,  V,  w 

AAA 

U,  V,  w 

u,  v,  W 

z 

“kn 

“kn 

6C 

Apa 

APj 

AP1/2 

At 

e 

9 

v 

vf 

P 


0) 


ft(t) 


flow  velocities/(a  fy,)  in  radial,  azinuthal,  and  axial 
directions,  respectively 

flow  perturbation  velocities/(a  ftQ)  in  r,  e,  and  z 
directions,  respectively 

radial  variation  of  u,  v,  w  ((la)  -  ( 1c ) ) 

unperturbed  flow  velocities/(a  fl  )  in  r,  e,  and  z  directions, 
respectively 

axial  coordinate  (z  =  0  at  base  of  cylinder) 

decay  rate/fiQ  of  k,n  mode  determined  by  numerical  simulation 

damping  factor  used  in  preprocessing  data  ((7)) 

complex  boundary  layer  thickness  occurring  in  the  viscous 
endwall  correction  to  the  eigenvalue  solution 

~  0 • b  Ap  ^  “  Ap  2/2 


=  p(0,c,t)  -  p(0,0,t) 

=  p(0,c/2,t)  -  p(0,0,t) 
time  interval 

amplitude  of  perturbation  to  angular  velocity  of  cylinder 

azimuthal  cylindrical  coordinate 

kinematic  viscosity  of  liquid 

frequency  in  Fourier  transform  plots  [Hz] 

density  of  liquid 

stream  function  in  Navier-Stokes  spin-up  flow  calculation  ((A. 2)) 

frequency/^,  of  perturbation  to  cylinder  spin 

value  of  to  at  which  peak  occurs  on  response  curve  in  Figure  3 
unperturbed  angular  velocity  of  cylinder 

forcing  function  of  perturbation  to  cylinder  spin 
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this  report  to  raise  specific  questions  or  discuss  the  topic, 
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